      SUBROUTINE CSTAB(A,NA,B,NB,F,NF,IOP,SCLE,DUMMY)
      IMPLICIT REAL*8(A-H,O-Z)
C
C
      DIMENSION A(1),B(1),F(1),DUMMY(1)
      DIMENSION NA(2),NB(2),NF(2),IOP(3),NDUM(2)
      DIMENSION IOPT(2)
      LOGICAL SYM
      COMMON/TOL/EPSAM,EPSBM,IACM
C
      N = NA(1)**2
      N1=N+1
C
      IF(IOP(2) .EQ. 0 ) GO TO 100
      CALL EQUATE(A,NA,DUMMY,NA)
      N2=N1+NA(1)
      N3=N2+NA(1)
      ISV =0
      ILV =0
      CALL EIGEN(NA(1),NA(1),DUMMY,DUMMY(N1),DUMMY(N2),ISV,ILV,V,DUMMY(N
     13),IERR)
C
      M=NA(1)
      IF(IERR .EQ. 0) GO TO 50
      CALL LNCNT(3)
      PRINT 25,IERR
   25 FORMAT(//' IN CSTAB, THE SUBROUTINE EIGEN FAILED TO DETERMINE THE
     1',I4,' EIGENVALUE FOR THE MATRIX A  AFTER 30 ITERATIONS')
      IERR=1
      CALL NORMS(M,M,M,A,IERR,BETA)
      BETA=2.*BETA
      GO TO 200
   50 CONTINUE
C
      BETA = 0.0
      DO 75 I = 1,M
      J = N1 + I - 1
      BETA1 = ABS(DUMMY(J))
      IF(BETA1 .GT. BETA)  BETA = BETA1
   75 CONTINUE
      BETA = SCLE*(BETA + .001)
      GO TO 200
C
  100 CONTINUE
      BETA = SCLE
  200 CONTINUE
C
      CALL TRANP(B,NB,DUMMY,NDUM)
      CALL MULT(B,NB,DUMMY,NDUM,DUMMY(N1),NA)
      CALL SCALE(DUMMY(N1),NA,DUMMY,NA,-2.0)
      CALL SCALE(A,NA,DUMMY(N1),NA,-1.0)
      J = -NA(1)
      NAX = NA(1)
      DO 225 I=1,NAX
      J = J+NAX+1
      K = N1+J-1
      DUMMY(K)=DUMMY(K)-BETA
  225 CONTINUE
      N2 = N1 + N
      SYM = .TRUE.
      IOPT(1)=0
C
      IF( IOP(3) .NE. 0 )  GO TO 300
      EPSA=EPSAM
      CALL BARSTW(DUMMY(N1),NA,A,NA,DUMMY,NA,IOPT,SYM,EPSA,EPSA,DUMMY(N2
     1))
      GO TO 350
  300 CONTINUE
      IOPT(2) = 1
      CALL BILIN(DUMMY(N1),NA,A,NA,DUMMY,NA,IOPT,ASCLE,SYM,DUMMY(N2))
  350 CONTINUE
C
      CALL EQUATE(B,NB,DUMMY(N1),NB)
      IOPT (1)= 3
      IAC =IACM
      N3 = N2 + NA(1)
      CALL SNVDEC(IOPT,NA(1),NA(1),NA(1),NA(1),DUMMY,NB(2),DUMMY(N1),IAC
     1,ZTEST,DUMMY(N2),DUMMY(N3),IRANK,APLUS,IERR)
      IF(IERR .EQ. 0 ) GO TO 400
      CALL LNCNT(5)
      IF(IERR .GT. 0 ) PRINT 360,IERR
      IF(IERR .EQ. -1) PRINT 370,ZTEST,IRANK
  360 FORMAT(//' IN CSTAB, SNVDEC HAS FAILED TO CONVERGE TO THE ',I4,' S
     1INGULARVALUE AFTER 30 ITERATIONS'//)
  370 FORMAT(//' IN CSTAB, THE MATRIX SUBMITTED TO SNVDEC USING ZTEST =
     1',E16.8,' IS CLOSE TO A  MATRIX OF LOWER  RANK '/' IF THE ACCURACY
     2 IAC IS REDUCED THE RANK MAY ALSO BE REDUCED'/' CURRENT RANK =',I4
     2)
      IF( IERR .GT. 0 ) RETURN
      NDUM(1) = NA(1)
      NDUM(2) =1
      CALL PRNT(DUMMY(N2),NDUM,'SGVL',1)
  400 CONTINUE
C
      CALL TRANP(DUMMY(N1),NB,F,NF)
      IF ( IOP(1) .EQ. 0 )   RETURN
      CALL LNCNT(4)
      PRINT 500
  500 FORMAT(//' COMPUTATION OF F MATRIX SUCH THAT A-BF IS ASYMPTOTICALL
     1Y STABLE IN THE CONTINUOUS SENSE '/)
      CALL PRNT(A,NA,' A  ',1)
      CALL LNCNT(4)
      PRINT 550,BETA
  550 FORMAT(//' BETA = ',E16.8/)
      CALL PRNT(B,NB,' B  ',1)
      CALL PRNT(F,NF,' F  ',1)
      CALL MULT(B,NB,F,NF,DUMMY,NA)
      CALL SUBT(A,NA,DUMMY,NA,DUMMY,NA)
      CALL PRNT(DUMMY,NA,'A-BF',1)
      N2 = N1+NA(1)
      N3 = N2+NA(1)
      ISV = 0
      ILV = 0
      CALL EIGEN(NA(1),NA(1),DUMMY,DUMMY(N1),DUMMY(N2),ISV,ILV,V,DUMMY(N
     13),IERR)
      M = NA(1)
      IF( IERR .EQ. 0 ) GO TO 600
      M = NA(1)-IERR
      CALL LNCNT(3)
      PRINT 25,IERR
  600 CONTINUE
      CALL LNCNT(4)
      PRINT 650
  650 FORMAT(//' EIGENVALUES OF A-BF'/)
  675 FORMAT(10X,2E16.8)
      CALL LNCNT(M)
      DO 700 I=1,M
      J = N1+I-1
      K = N2+I-1
      PRINT 675,DUMMY(J),DUMMY(K)
  700 CONTINUE
C
      RETURN
      END
